Homework 7.2: Dynamics of a Chemostat

a) d$c$/d$t$

b) d$n$/d$t$

c) nondimensionalization

\begin{align} \cfrac{\mathrm{d}c}{\mathrm{d}t} &= \beta \cfrac{n/k}{1+n/k} c - \cfrac{Q}{V}c \\[0.5em] \cfrac{\mathrm{d}n}{\mathrm{d}t} &= \cfrac{Qn_0}{V} - \cfrac{Qn}{V} - \nu \beta \cfrac{n/k}{1+n/k} c \end{align}

$\hspace{15em}\mathrm{Nondimensionalizing... }\\[0.5em]$ \begin{align} \cfrac{c_d}{t_d}\cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}t} &= \beta \cfrac{n_d\tilde{n}/k}{1+n_d\tilde{n}/k} c_d \tilde{c} - \cfrac{Q}{V}c_d \tilde{c} \\[0.5em] \cfrac{n_d}{t_d}\cfrac{\mathrm{d}\tilde{n}}{\mathrm{d}t} &= \cfrac{Qn_0}{V} - \cfrac{Qn_d\tilde{n}}{V} - \nu \beta \cfrac{n_d\tilde{n}/k}{1+n_d\tilde{n}/k} c_d \tilde{c} \\[2.5em] \end{align}

$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{n_d = k, \hspace{0.5em} t_d = \cfrac{V}{Q}} \\[1.0em] \end{align}

\begin{align} \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}t} &= \beta \cfrac{V}{Q} \cfrac{\tilde{n}}{1+\tilde{n}} \tilde{c} - \tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{n}}{\mathrm{d}t} &= \cfrac{n_o}{k} - \tilde{n} - \nu \beta \cfrac{V}{Q} \cfrac{c_d}{k} \cfrac{\tilde{n}}{1+\tilde{n}}\tilde{c} \\[0.5em] \end{align}

$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{ \tilde{\beta} = \beta \cfrac{V}{Q}, \hspace{0.8em} \tilde{n_0} = \cfrac{n_0}{k}, \hspace{0.8em} c_d = k } \end{align}

\begin{align} \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}t} &= \tilde{\beta} \cfrac{\tilde{n}}{1+\tilde{n}} \tilde{c} - \tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{n}}{\mathrm{d}t} &= \tilde{n_0} - \tilde{n} - \nu \tilde{\beta} \cfrac{\tilde{n}}{1+\tilde{n}}\tilde{c} \end{align}

$\hspace{15em}\mathrm{Rewriting... }$ \begin{align} \boxed{ \cfrac{\mathrm{d}c}{\mathrm{d}t} = \beta \cfrac{n}{1+n} c - c \\[0.5em] \cfrac{\mathrm{d}n}{\mathrm{d}t} = n_0 - n - \nu \beta \cfrac{n}{1+n}c } \end{align}

d) analysis

$\hspace{10em}\text{We find the nullclines and fixed points by setting the derivatives to zero:}$ \begin{align} \mathrm{nullcline}_1. \hspace{0.5em} &c = 0, && n = n_0 \\[0.5em] \mathrm{nullcline}_2. \hspace{0.5em} &c = \cfrac{(n_0 - n)(1+n)}{\nu \beta n}, && n = \cfrac{1}{\beta-1} \\[0.5em] \mathrm{fp}_1. \hspace{0.5em} &c_{st} = 0, &&n_{st} = n_0 \\[0.5em] \mathrm{fp}_2. \hspace{0.5em} &c_{st} = \cfrac{1}{\nu}\left(n_0 - \cfrac{1}{\beta-1} \right), && n_{st} = \cfrac{1}{\beta - 1} \\[0.5em] \end{align}

$\hspace{16em}\text{Note that for $\mathrm{fp}_2$ to be valid, $\boxed{\beta > 1, \hspace{0.5em} n_0 \geq \cfrac{1}{\beta-1}}$} $

$\hspace{10em}\text{We can then write out the linear stability matrix:}$

\begin{align} A = \begin{pmatrix} \beta \cfrac{n}{1+n} - 1 & \beta c \cfrac{1}{(1+n)^2} \\ -\nu \beta \cfrac{n}{1+n} & -1 -\nu \beta c \cfrac{1}{(1+n)^2} \end{pmatrix} \end{align}

$\hspace{10em}\text{Using Mathematica, the eigenvalues are:}$

\begin{align} \lambda_i = -1, \hspace{2em} \lambda_{ii} = \cfrac{(\beta-1)n^2 + (\beta-2)n - (\beta c \nu + 1) }{(1+n)^2} \end{align}

At the fixed point $\mathrm{fp}_1$ \begin{align} \lambda_{ii} &= \cfrac{n_0 (\beta-1) - 1}{1+n_0} < 0 \\[0.5em] &\text{when either 1. } \beta < 1,\text{or 2. } n_0 < \cfrac{1}{\beta-1} \hspace{0.1em}\mathrm{and}\hspace{0.1em} \beta > 1. \end{align}

At the fixed point $\mathrm{fp}_w$ \begin{align} \lambda_{ii} &= 1-(\beta-1)n_0 < 0 \\[0.5em] &\mathrm{if}\hspace{0.6em} n_0 > \cfrac{1}{\beta - 1}\hspace{12em} \end{align}

Note that for $\mathrm{fp}_2$ we take $\beta$ to be greater than 1 since its existence is predicated on this condition.
We can break up our analysis into two cases:

$\beta \leq 1: $

Production is not fast enough relative to how quickly $c$ is being flushed out. Only one fixed point exists where the bacterial concentration decays to zero, and the nutrient concentration steadies to $n_0$. This fixed point has two negative eigenvalues (look at numerator of $\lambda_{ii}$), and is thus stable.

$\beta > 1: $

When production is fast enough, we also need the influx nutrient concentration to be sufficiently high, $n_0 \geq 1/(\beta-1)$ to reach a nonzero bacterial concentration. When both conditions are satisfied, fixed point 2 (nonzero $c$) is stable, and fixed point 1 (zero $c$) is unstable (condition 2 is not satisfied since the inequality is flipped).
When $n_0$ is not sufficiently high, fixed point 1 (zero $c$) is again stable (condition 2 is satisfied once more).

First some colors.

Define some functions next.

Nullclines

Let's call all our functions! I am setting $\beta > 1$ and $n_0 > 1/(\beta-1)$, and we see exactly what we expect, with the entire positive domains swirling into our steady state.

Let's look at the three cases outlined above. On the left, we have the case when both conditions are satisfied. In the middle, n0 is not sufficiently large. On the right, $\beta$ is less than 1. In both cases we venture into the negatives.

Note the subtle changes in the stabilities. I plot things in the negatives to get a better idea of what's happening with the math. Big red lines are drawn on the zeroes.

Trajectories

Let's plot some trajectories. Note that my initial conditions are thrown in randomly, but the random distribution is over a log scale to get better resolution at lower concentrations.

Here's a quick snapshot of the different behaviors.

Observations

Toggling $\beta$ and $n_0$, we get precisely the trend we expect, where $c$ falls in the regime of low $\beta < 1$, and rises to a steady state when both n0 and β are large. Increasing nu decreases our steady state $c$, and increasing n0 increases our steady state $c$. Since nu is harder to tune in real life, what we have to do to experimentally is:

The above inequalities yield the following. \begin{align} Q & < \beta V \\[0.5em] Q & < \beta V \cfrac{1}{1 + k/n_0} \end{align} Note that the second condition yields the stricter inequality. Thus the chosen inlet flow rate should be low enough to give the system enough time to build up sufficient quantities of $c$ before it gets flushed out. $$\boxed{Q < \beta V \cfrac{1}{1 + k/n_0}}$$

Experimental design

Since $\beta$, $k$, and $\nu$ are intrinsic to the nutrients and bacteria, we can tune $n_0$ to give the desired steady-state $c$.

Note: we would have to solve Q and n_0 simultaneously. The general idea in practice is to first try a pair of ($n_0$, Q) with small Q/V relative to $\beta$, and $n_0$ close to $k$ to then yield a desirable $\tilde{c_{st}}$, then adjust Q according to the final inequality. For some experimental values of $V, k, \beta, \nu$, we can also computationally find the perfect pair.

The color represents the [c]. The black line represents $Q = \beta V / (1+k/n_0)$. All points below this curve satisfy both relations. We see that above the curve, [c] dips into the negatives. We thus see that we can pass in the fixed values of all the other parameters, locate the desired $c$ on the color-bar, then pick an ($n_0$, $Q$) pair. $Q$ will partly control the timescales of the experiment (how long it takes to reach this steady state), so the fact we have a log-scale to choose from is very nice.

Concluding thoughts

In this problem, I really enjoyed how elegant the results were even though it was a relatively simple system, because I felt like I understood why every parameter would be important in the context of an experiment. Chemostat, thank you for the chat! And reader, as always, thank you for reading!